#include <string.h>
#include <cstring>
#include <time.h>
#include <math.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <sstream>


#include <omp.h>

#ifdef _WIN32
#include <direct.h>
#define vbmkdir(x) _mkdir(x)
#else
#define vbmkdir(x) mkdir(x,S_IRWXU|S_IRGRP|S_IXGRP )
#endif 

#include "DDP.h"
#include "HapAssembler.h"

using namespace std;

////////////////// for initialization ////////////////////
int		HAP_BL_LENGTH	= 7;
int		HAP_BL_START	= 0;
int		HAP_BL_END		= 10000;
int		HAP_BL_MAX	= 10000;
int		HAP_BL_MIN	= 0;
int		HAP_NUM_BURN_IN = -1; //5000; //50000;
int		HAP_NUM_CUM = -1; // 3000; // 50000;
int		HAP_THINING = -1;
int		HAP_CONPARAM = -1;
int		HAP_CONVCHECK = -1;
const char *HAP_INI_FILE;
char *HAP_OUTDIR = 0;
const char *HAP_MA = NULL;
void copyAncestor(DDP Haplo, int blockNum, FILE * myfile);
int deleteFile(char *filename);


///////////////////////////////////////////
/// parse argument
///////////////////////////////////////////
void parseArg( int argc, char *argv[] )
{
	if ( argc < 2 )
	{
		printf("[usage]: Haplo.exe inputfile\n");
		exit(-1);
	}

	HAP_INI_FILE = argv[1];
	int counter = 1;

	while (++counter < argc) 
	{

	    if(strcmp(argv[counter], "-ma") == 0){//is MA is given as Input
		counter++;
		HAP_MA = argv[counter];
	    }
	    else if ( strcmp( argv[counter], "-bl" ) == 0 ) // block length in Partition-Ligation for long sequences
		{
			counter++;
			HAP_BL_LENGTH = atoi( argv[counter] );			
		} 
		else if ( strcmp ( argv[counter], "-dir" ) == 0 )	// output dir name
        {
              counter++;
              HAP_OUTDIR = argv[counter] ;
        }
		else if ( strcmp ( argv[counter], "-cp" ) == 0 )	// flag for resampling DP hyper parameters
		{
			counter++;
			HAP_CONPARAM = atoi( argv[counter] );			
		}
		else if ( strcmp ( argv[counter], "-conv" ) == 0 )	// flag for checking convergence or doing fixed # of iter
		{
			counter++;
			HAP_CONVCHECK = atoi( argv[counter] );			
		}

		else if ( strcmp ( argv[counter], "-nbr" ) == 0 )	// number of burnin iterations
		{
			counter++;
			HAP_NUM_BURN_IN = atoi( argv[counter] );			
		}
		else if ( strcmp ( argv[counter], "-nc" ) == 0 )	// number of cumulative iterations
		{
			counter++;
			HAP_NUM_CUM = atoi( argv[counter] );		
		}
		else if ( strcmp ( argv[counter], "-nthin" ) == 0 )	// thining interval
		{
			counter++;
			HAP_THINING = atoi( argv[counter] );		
		}
		else if ( strcmp ( argv[counter], "-mbl" ) == 0 )	// max SNP index
		{
			counter++;
			HAP_BL_MAX = atoi( argv[counter] );			
		}
		else if ( strcmp ( argv[counter], "-mnbl" ) == 0 )	// min SNP index
		{
			counter++;
			HAP_BL_MIN = atoi( argv[counter] );			
		}
		else 
		{
			printf( "Unknown option: [%s]\n", argv[counter]);
//		    exit(-1);
		}		
	} // while
}

///////////////////////////////////////////
/// Load ini file
///////////////////////////////////////////
vector<string> LoadIniFile( const char* inifile )
{
	char	tmp[500];
	int		nfile;
	vector<string>	filelist;
	
	FILE *fp = fopen( inifile, "r" );
	if ( fp == 0 )
	{
		printf( "File Not Exist: %s\n", inifile );
		exit(-1);
	}

	fscanf( fp, "%d", &nfile );			// line 1. number of files (populations)
	for (int i=0; i < nfile; i++)		// line 2-n. filename for each population data
	{
		fscanf( fp, "%s", tmp );
		filelist.push_back( tmp );
		printf( "%s\n", filelist[i].c_str());
	}
	
	fclose(fp);
	
	return filelist;
}

///////////////////////////////////////////
/// Make output directory
///////////////////////////////////////////
int MakeDir( const char *dir )
{
	string dname( dir );

	int sep, sti=0;
	string rem(dname);
	string sub;
	while (1)
	{
		sep = rem.find( "/" );
		if ( sep == -1 )
			return 1;

		sti += sep+1;
		sub = dname.substr(0, sti);
		vbmkdir( sub.c_str() );
		rem = rem.substr( sep+1, rem.length() );
	}

	return 1;
}

///////////////////////////////////////////
int main( int argc, char *argv[] )
{
int numthreads;
double wtime;
#ifdef _OPENMP
   numthreads=omp_get_max_threads();
   wtime = omp_get_wtime();
#else
   numthreads=1;
#endif

	/////////////////////////////////
	// parse arguments
	parseArg( argc, argv );

	///////////////////////////////////////////////////
	///////// parse file path ////////////
	string dbname( HAP_INI_FILE );

	int pos = dbname.rfind("/");
	if ( pos > 0 )
			dbname = dbname.substr( pos+1, dbname.length() );
	int dot = dbname.rfind( "." );
	dbname = dbname.substr(0, dot);



	////////////////////////////////////////////
	//creating output directory

printf("Number of Threads is %d\n", numthreads);
	ostringstream buffer;
        buffer << wtime;
	string str = buffer.str();
	char dir[150] = "./output/";
	//strcat(dir,str.c_str());
	//strcat(dir,".");
	strcat(dir,dbname.c_str());
	strcat(dir,"/");

	////////////////////////////////////////////
	// creating m_A path
	char path[150];
	strcpy(path, dir);
	strcat(path,"m_A");
	char * ANCESTOR_OUT = path;


	////////////////////////////////////////////
	string outdir( dir );
	if ( HAP_OUTDIR != 0 )
		outdir = HAP_OUTDIR;

	MakeDir( outdir.c_str() );


	/////////////////////////////////
	// Read ini file
	//vector<string>	filelist = LoadIniFile( HAP_INI_FILE );
//govind - changed file entering to single file
        vector<string> filelist;
	filelist.push_back(HAP_INI_FILE);

	//////////////////////////////////
	// Load DB
	GenoHaploDB	DB;
	DB.m_nMaxT = HAP_BL_MAX;
	DB.m_nMinT = HAP_BL_MIN;

        bool isM_A = false;//whether M_A file is specified or not
	if(HAP_MA != NULL){
	    isM_A = true;
	    filelist.push_back(HAP_MA);
	    DB.LoadMA(filelist, DB.m_A);
	}

	DB.LoadData( filelist, isM_A );
	HAP_BL_MAX = DB.m_nMaxT;
	int		I = DB.m_numTotalI;
	int		numT = DB.m_numT;
	int		nGroups = filelist.size();

	//////////////////////////////////
	clock_t start, finish;
	start = clock();




	////////////////////////////////////////////////
	////	create block-assembler to ligate blockwise result
	HapAssembler Assembler( outdir, dbname ) ;

	vector<haplo2_t> &arrData = Assembler.m_arrData;
	vector<haplo2_t> &arrPred = Assembler.m_arrPairwisePred;


	int res = numT % HAP_BL_LENGTH;
	int	numBlocks = (numT - res) / HAP_BL_LENGTH;
	if ( res !=0 ) 
			numBlocks++;

	
	int		nstart, nend, nlength;

	FILE * myfile;
	bool copyMA = false; //whether we want to copy MA or not
	if(copyMA){	
	    deleteFile(ANCESTOR_OUT);
	    myfile = fopen(ANCESTOR_OUT, "w");
	}
	DDP Haplo;

	////////////////////////////////////////////////
	////	create HDP instance to phase each block
	

	///////////////////////////////////////
	///		1. Block-wise inference  using HDP haplotyper (Partition)
	///////////////////////////////////////
#pragma omp parallel for  shared(DB) private(nend, nstart, nlength) lastprivate(Haplo)
	for ( int nb = 0; nb < numBlocks; nb++ )
	{
	    Haplo.Init( &DB, 0, numT );
	    //Burning is transitory starting time of markov chain
	    if ( HAP_NUM_BURN_IN > 0 )
		Haplo.m_nBurninIteration = HAP_NUM_BURN_IN;
	    if ( HAP_NUM_CUM > 0 )
		Haplo.m_nCumIteration = HAP_NUM_CUM;
	    if ( HAP_THINING > 0 )
		Haplo.m_nThining = HAP_THINING;
	    if (  HAP_CONPARAM > 0 )
		Haplo.m_doConparam = HAP_CONPARAM;
	    if (  HAP_CONVCHECK > 0 )
		Haplo.m_bCheckConvg = HAP_CONVCHECK;

	    Assembler.Initialize( &DB, &Haplo );

	    nend = MIN( ( nb + 1 ) * HAP_BL_LENGTH, numT );
	    nstart = MAX( nend - HAP_BL_LENGTH, 0);
	    nlength = nend - nstart;
            vector<vector<unsigned char > > lma;
	    Haplo.m_A = DB.m_A[nb];
	    //lma = DB.m_A[nb];
	    if ( HAP_BL_LENGTH != nlength )
	    {	
		HAP_BL_LENGTH = nlength;
		cout << "block length not equal to " << HAP_BL_LENGTH << endl;
	    }	

	    Haplo.InferHaplotypes( &DB, nstart, nend );

	    // save intermediate result ( just for checking, can be commented out )
	    Haplo.Save( dbname.c_str(), outdir.c_str(), nstart, nend );

	    // add the blockwise result to the output array
	    if(copyMA){
		copyAncestor(Haplo, nb, myfile);
	    }
	}
	fclose(myfile);

	for ( int nb = 0; nb < numBlocks; nb++ )
	{
		nend = MIN( ( nb + 1 ) * HAP_BL_LENGTH, numT );
		nstart = MAX( nend - HAP_BL_LENGTH, 0);
		nlength = nend - nstart;
		haplo2_t	*phh = new haplo2_t;
		DB.AllocHaplo2( phh, I, nlength );
		DB.CopyHaplo( nstart, nend, phh );
		arrData.push_back( *phh );
        }
	///////////////////////////////////////////////
	///////// 2. Ligation 
	///////////////////////////////////////////////
	vector<int> T1( numBlocks, HAP_BL_LENGTH );		// length of each block
	vector<int> nOverlap( numBlocks, 0 );			// overlap between adjacent block
	if ( res != 0 )
		nOverlap[ numBlocks - 2 ] = HAP_BL_LENGTH - res;

	Assembler.LigateHaplotypes( arrData, T1, nOverlap, arrPred );


	//////////// save final result ///////////////
	// output filename
	string ofilename( outdir );
	ofilename += dbname + "_final.txt";
	Assembler.SaveResult( ofilename.c_str() );
	
	finish = clock();
	
//	printf("Elapsed time: %.2f min\n", (double)(finish-start)/CLOCKS_PER_SEC/60.0 );
	printf("Elapsed time: %.2f min\n", (double)(finish-start)/1000/60.0 );

#ifdef _OPENMP
	wtime = omp_get_wtime()- wtime;
	fprintf (stderr, "  Elapsed wall clock time (seconds) %f\n", wtime );
	fflush(stderr);
#endif
	return 1;
}

void copyAncestor(DDP Haplo, int blockNum, FILE * myfile){
    int totalAncestor = Haplo.m_A.size();
    int nstart = Haplo.m_nBlockStart;
    int nend = Haplo.m_nBlockEnd;
    ostringstream myStream;
    char buffer[10000];

    sprintf(buffer, "%d:",blockNum );

    for(int k=0; k < totalAncestor; k++){
	    for(int tt=0; tt < nend - nstart; tt++){
		    sprintf(buffer + strlen(buffer),"%u", Haplo.m_A[k][tt]);
	    }
            sprintf(buffer + strlen(buffer), ":");
    }

    sprintf(buffer + strlen(buffer), "\n");
    fprintf(myfile, "%s", buffer);
}

int deleteFile(char *filename){
    ifstream fin(filename);
    if (fin)
    {
        char cmd[50];
        remove(filename);
        fin.close();
    }
    else
    {
        cout << "Could not find file!" << endl;
    }
    return 0;
}

